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ABSTRACT 

Efficient algorithms are known for many operations on trun- 
cated power series (multiplication, powering, exponential, 
...). Composition is a more complex task. We isolate a 
large class of power series for which composition can be 
performed efficiently. We deduce fast algorithms for con- 
verting polynomials between various bases, including Euler, 
Bernoulli, Fibonacci, and the orthogonal Laguerre, Hermite, 
Jacobi, Krawtchouk, Meixner and Meixner-PoUaczek. 

Categories and Subject Descriptors: 
1.1.2 [Computing Methodologies]: Symbolic and Alge- 
braic Manipulation - Algebraic Algorithms 

General Terms: Algorithms, Theory 

Keywords: Fast algorithms, transposed algorithms, basis 
conversion, orthogonal polynomials. 

1. INTRODUCTION 

Through the Fast Fourier Transform, fast polynomial mul- 
tiplication has been the key to devising efflcient algorithms 
for polynomials and power series. Using techniques such 
as Newton iteration or divide-and-conquer, many problems 
have received satisfactory solutions: polynomial evaluation 
and interpolation, power series exponentiation, logarithm, 
. . . can be performed in quasi-linear time. 

In this article, we discuss two questions for which such fast 
algorithms are not known: power series composition and 
change of basis for polynomials. We isolate special cases, 
including most common families of orthogonal polynomials, 
for which our algorithms reach quasi-optimal complexity. 

Composition. Given a power series g with coefficients in 
a field K, we first consider the map of evaluation at g 



Eva\m,n{-,g) ■■ Ae 



A{g) mod x" £ 



Here, K[a;]m is the m-dimensional K-vector space of polyno- 
mials of degree less than m. We note Eval„ for Eval„,„. 

To study this problem, as usual, we denote by M a mul- 
tiplication time function, such that polynomials of degree 
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less than n can be multiplied in M(n) operations in K. We 
impose the usual super-linearity conditions of [17, Chap. 8]. 
Using Fast Fourier Transform algorithms, M(n) can be taken 
in 0(nlog(n)) over fields with suitable roots of unity, and 
0(n log(7i) log log(n)) in general [31, 14]. 

If g{0) = 0, the best known algorithm, due to Brent and 
Kung, uses 0{y^nlogn M(n)) operations in K [11]; in small 
characteristic, a quasi-linear algorithm is known [5]. There 
are however special cases of power series g with faster algo- 
rithms: evaluation ai g = Xx takes linear time; evaluation 



at g 



requires no arithmetic operation. A non-trivial 



example is g = x + a, which takes time 0(M(n)) when 
the base field has characteristic zero or large enough [1]. 
Brent and Kung [11] also showed how to obtain a cost in 
0(M(n) log(n)) when g is a polynomial; this was extended 
by van der Hoeven [22] to the case where g is algebraic over 
K(a;). In §2, we prove that evaluation at <; = exp(a::) — 1 and 
at g = log(l -I- x) can also be performed in 0{M{n) log(n)) 
operations over fields of characteristic zero or larger than n. 
Using associativity of composition and the linearity of 



the map Evain 



we show in §2 how to use these spe- 



cial cases as building blocks, to obtain fast evaluation al- 
gorithms for a large class of power series. This idea was 
first used by Pan [28], who applied it to functions of the 
form {ax + b)/{cx + d). Our extensions cover further exam- 
ples such as 2x/(l + x)^ or (1 — \/l — x'^)/x, for which we 
improve the previously known costs. 

Bivariate problems. Our results on the cost of evalu- 
ation (and of the transposed operation) are applied in §3 
to special cases of a more general composition, reminiscent 
of umbral operations [30]. Given a bivariate power series 
F = "^2 x} (.ii^W , we consider the linear map 

Eval„(.,F,i) : (ao, . . . ,a„_i) >-» ^Ci(a;)«j mod a;". 

For instance, with 

1 



F = 



1 - tg{x) 



Y,9{^yt' 



J>0 



this is the map Eva In (.,3) seen before. For general F, the 
conversion takes quadratic time (one needs n^ coefficients 
for F). Hence, better algorithms can only been found for 
structured cases; in §3, we isolate a large family of bivari- 
ate series F for which we can provide such fast algorithms. 
This approach follows Frumkin's [16], which was specific to 
Legendre polynomials. 

Change of basis. Our framework captures in particular 
the generating series of many classical polynomial families. 



for which it yields at once conversion algorithms between 
the monomial and polynomial bases, in both directions. 

Thus, we obtain in §4 change of basis algorithms of cost 
only 0(M(n)) for all of Jacobi, Laguerre and Hermite or- 
thogonal polynomials, as well as Euler, Bernoulli, and Mott 
polynomials (see Table 3). These algorithms are derived in 
a uniform manner from our composition algorithms; they 
improve upon the existing results, of cost 0(M(n) log(n)) or 
0(M(n) log^(n)) at best (see below for historical comments). 

We also obtain 0(M(n) log(n)) conversion algorithms for a 
large class of Sheffer sequences [30, Chap. 2], including actu- 
arial polynomials, Poisson-Charlier polynomials and Meix- 
ner polynomials (see Table 4). 

Transposition. A key aspect of our results is their heavy 
use of transposed algorithms. Introduced under this name 
by Kaltofen and Shoup, the transposition principle is an 
algorithmic theorem with the following content: given an 
algorithm that performs an r x s matrix- vector product b i-^ 
Mb, one can deduce an algorithm with the same complexity, 
up to 0{r + s) operations, and that performs the transposed 
matrix-vector product c t-^ M^c. In other words, this relates 
the cost of computing a K-linear map / : V — > W^ to that 
of computing the transposed map /* : W* -^ V* . 

For the transposition principle to apply, some restrictions 
must be imposed on the computational model: we require 
that only linear operations in the coefficients of b are per- 
formed (all our algorithms satisfy this assumption). See [12] 
for a precise statement, Kaltofen's "open problem" [23] for 
further comments and [7] for a systematic review of some 
classical algorithms from this viewpoint. 

To make the design of transposed algorithms transparent, 
we choose as much as possible to describe our algorithms 
in a "functional" manner. Most of our questions boil down 
to computing linear maps K[a;]m -^ ]K[a;]„; expressing algo- 
rithms as a factorization of these maps into simpler ones 
makes their transposition straightforward. In particular, 
this leads us to systematically indicate the dimensions of 
the source (and often target) space as a subscript. 

Previous work. The question of efficient change of ba- 
sis has naturally attracted a lot of attention, so that fast 
algorithms are already known in many cases. 

Gerhard [18] provides 0(M(n) log(n)) conversion algori- 
thms between the falling factorial basis and the monomial 
basis: we recover this as a special case. The general case of 
Newton interpolation is discussed in [6, p. 67] and developed 
in [9]. The algorithms have cost 0(M(n) log(n)) as well. 

More generally, if [Pi) is a sequence of polynomials satis- 
fying a recurrence relation of fixed order (such as an orthog- 
onal family), the conversion from {Pi) to the monomial basis 
(a;') can also be computed in 0(M(n) log(n)) operations: an 
algorithm is given in [29], and an algorithm for the trans- 
pose problem is in [15]. Both operate on real or complex 
arguments, but the ideas extend to more general situations. 
Alternative algorithms, based on structured matrices tech- 
niques, are given in [21]. They perform conversions in both 
directions in cost 0(M(n) log^(n)). 

The overlap with our results is only partial: not all fam- 
ilies satisfying a fixed order recurrence relation fit in our 
framework; conversely, our method applies to families which 
do not necessarily satisfy such recurrences (the work-in- 
progress [8] specifically addresses conversion algorithms for 
orthogonal polynomials). 

Besides, special algorithms are known for converting be- 



tween particular families, such as Chebyshev, Legendre and 
Bezier [27, 4], with however a quadratic cost. Floating-point 
algorithms are known as well, of cost 0{n) for conversion 
from Legendre to Chebyshev bases [2] and 0(nlog(n)) for 
conversions between Gegenbauer bases [25], but the results 
are approximate. Approximate conversions for the Hermite 
basis are discussed in [26], with cost 0(M(n) log(n)). 

Note on the base field. For the sake of simplicity, in all 
that follows, the base field is supposed to have character- 
istic 0. All results actually hold more generally, for fields 
whose characteristic is sufficiently large with respect to the 
target precision of the computation. However, completely 
explicit estimates would make our statements cumbersome. 

2. COMPOSITION 

Associativity of composition can be read both ways: in 
the identity A[f o (?) = ^(/) ° 9, f is either composed on 
the left of g or on the right of A. In this section, we discuss 
the consequences of this remark. We first isolate a class of 
operators / for which both left and right composition can 
be computed fast. Most results are known; we introduce 
two new ones, regarding exponentials and logarithms. Us- 
ing these as building blocks, we then define composition se- 
quences, which enable us to obtain more complex functions 
by iterated compositions. We finally discuss the cost of the 
map Evain and of its inverse for such functions, showing how 
to reduce it to 0{M{n)) or 0(M(n) log(n)). 

2.1 Basic Subroutines 

We now describe a few basic subroutines that are the 
building blocks in the rest of this article. 

Left operations on power series. In Table 1, we list 
basic composition operators, defined on various subsets of 
]K[[a;]]. Explicitly, any such operator o is defined on a domain 
dom(o), given in the third column. Its action on a power 
series g £ dom(o) is given in the second column, and the 
cost of computing o{g) mod a;" is given in the last column. 



Operator Action 

Aa (add) 
Ma (mul) 
Pfe (power) 

Rk,a,r (root) 

Inv (inverse) 1/g 
E (exp.) exp(g) 



Domain 



Cost 



a + g 



K[[x]] 1 

K[[x]] n 

K[[x]] 0(logfc + M(n)) 
a''x''''il + xK[[x]]) 0(M(n)) 

]K*-ha:]K[[a:]] 0(M(n)) 

1 xK[[x]] 0(M(n)) 

L(log.) log(l+ff) xK[[x]] 0(M(n)) 

Table 1: Basic Operations on Pow^er Series 

Some comments are in order. For addition and multipli- 
cation, we take a G K and A in K*. To lift indeterminacies, 
the value of Rk.a,r{g) is defined as the unique power series 
with leading term ax'' whose fcth power is g; observe that to 
compute Rk,a,r{g) mod x", we need g modulo x"'^'~^''~^' as 
input. Finally, we choose to subtract 1 to the exponential so 
as to make it the inverse of the logarithm. All complexity re- 
sults arc known; they are obtained by Newton iteration [10]. 

Right operations on polynomials. In Table 2, we de- 
scribe a few basic linear maps on K[x]m (observe that the di- 
mension m of the source is mentioned as a subscript). Their 
action on a polynomial 



A{x) = ao H + a„ 



-ix""-^ e K[a;l. 



Name 



Notation 



Action 



Cost 



Powering Powers, fe A{x ) 

Reversal Rev™' x"'-'^A{l/x) 

Mod modm,n A mod x" 

Scale ScaleA,m A{\x) 0{m) 

Diagonal Am{-,Si) "^aiSiX^ m 

Multiply Mulm,n(.,P) ^Pmoda;" M(max(n,m)) 

Shift Shiftg,™ A{x + a) M(m) + 0(m) 

Table 2: Basic Operations on Polynomials 

is described in the third column. In the case of powering, 
it is assumed that k £ N>o- Here and in what follows, we 
freely identify K[a;]m and EC", through the isomorphism 



^a, 



x' G ¥.\xl 



(ao, • • • ,am-i) G K" 



All of the cost estimates are straightforward, except for the 
shift, which, in characteristic 0, can be deduced from the 
other ones by the factorization [1]: 

Shifta.m = Am(Rev„(MuU,,„(Rev™(Am(., i!)), P)), 1/j!), 

where P is the polynomial X]r=o ^'a^V*'- We continue with 
some equally simple operators, whose description however 
requires some more detail. For k G N>o, any polynomial A 
in 'K[x] can be uniquely written as 

A{x) = Ao/k{x'') + Aiik{x'')x + ■■■+ Ak-x/k{x'')x'''^. 

Inspecting degrees, one sees that if A is in K[a;]m, then Ai/^ 
is in ]K[a;]mi, with 

I 1 if i < m mod fc, , , 

mi=[m/k\+l - . (1) 

I otherwise. 

This leads us to define the map Split^ j, : 

A G K[x\m >-> {Ao/k, ■■■, Ak-i/k) e K[x]m„ X • • • X K[x]m^_^ . 

It uses no arithmetic operation. We also use linear com- 
bination with polynomial coefficients. Given polynomials 
Go, . . . , Gfc_i in K[a::]„i, we denote by 

Combm(.,Go, . . . ,Gfe_i) : K[x]^ — > K[x]m 

the map sending {Ao, . . . , Ak-i) G K[x]J^ to 

AoGo H h Ak-iGk-i mod x"" G ]K[a;]m. 

It can be computed in 0(fcM(m)) operations. Finally, we ex- 
tend our set of subroutines on polynomials with the following 
new results on the evaluation at exp(a;) — 1 and log(l + x). 

Proposition 1. The maps 

Exp^^ : A G K[x]m 1-^ ^(exp(a::) — 1) mod x" G K[a::]„, 

Log„_„ : A G K[a:;]„ i-^ ^(log(l + x)) mod x" G K[a;]„ 

can be computed in 0(M(n) log(n)) arithmetic operations. 

Proof. We start by truncating A modulo x" , since 

Exp„,„(A) = Exp„ „(^ mod x"). 

After shifting by —1, we are left with the question of evalu- 
ating a polynomial in K[a;]„ at X^^^^a;'/*!- Writing its ma- 
trix shows that this map factors as A„(MultiEvalJj(.), 1/i!), 
where MultiEval„ is the map 

A G K[2;]„ ^ (A(0), ...,A{n- 1)) G K". 



To summarize, we have obtained that 

Exp„,„(^) = A„(MultiEval^Shift_i,„(mod„,„(A))), 1/i!). 

Using fast transposed evaluation [13, 7], Exp^ n(^) can thus 
be computed in 0(M(n) log(n)) operations. Inverting these 
computations leads to the factorization 

Log„ „(A) = Shifti.„(lnterpJ,(A„(mod„,„(A),i!))), 

where lnterp„ is interpolation at 0, ... ,n — 1. Using algo- 
rithms for transpose interpolation [24, 7] , this operation can 
be done in time 0{M{n) log(n)). D 

2.2 Associativity Rules 

For each basic power series operation in Table 1, we now 
express E\/a\m,n{A,o{g)) in terms of simpler operations; we 
call these descriptions associativity rules. We write them 
in a formal manner: this formalism is the key to automati- 
cally design complex composition algorithms, and makes it 
straightforward to obtain transposed associativity rules, re- 
quired in the next section. Most of these rules are straight- 
forward; care has to be taken regarding truncation, though. 

Scaling, Shift and Powering. 

Eval„,„(^, Mx{g)) = EvaU,„(ScaleA,m(^),5), (Ai) 

Eva\^,„{A,Aa{g)) = EvaU,„(Shifta,™(A),g), (A2) 

E\/a\,n,n{A,Pk{g)) = E\/a\k(m-i)+i,n{Power„i^k{A),g). (A3) 

Inversion. From A(l/g) = (Re\/m{A)){g)/g"^^^ and writ- 
ing h — g^^"^ mod x", we get 

E\/a\m,n(A,\r\\/{g)) = Mul„,„(Evalm,„(ReVm(A),5), ft), 

(A4) 
Root taking. For g and h in K[[a;]], if 5 — h'', one has 
A{h) = Ao/k{g) + Ai/k{g)h + ■■■ + Ak-i/k{g)h''-\ We de- 
duce the following rule, where the indices nii are defined in 
Equation (1). 

hi = /i' mod a::" for < i < fe 

Ao,...,Ak-i^Sp\\t^^k{A) (As) 

Bi = E\/a\m. ^„{Ai,g) for < i < fc 

E\/a\m,n{A, Rk,a,r{g)) = Comb„(Bo,. . . ,Bk-i, 1, . . . ,hk-i). 

Exponential and Logarithm. 

EvaU.„(A, E{g)) = Eval„(Exp„_„(^), 5), (Ag) 

EvaU,„(A, L{g)) = Eval„(Log„^„(A),5). (A7) 

2.3 Composition sequences 

We now describe more complex evaluations schemes, ob- 
tained by composing the former basic ones. 

Definition 1. Let O be the set of actions from Table 1. 
A sequence o — (oi,...,o_l) with entries in O is defined 
at a series g G K[[a;]] tf g is in dom(oi), and for i < L, 
Oi-i(- • • Oi(p)) is in dom(oi). It is a composition sequence 
if it is defined at x; in this case, o computes the power series 
gi, . . . ,gL, with go = X and gi — Oi{gi-i); it outputs gL. 

Examples. As mentioned in [28], the rational scries g = 
{ax + b)/{cx + d) G ]K[[a::]], with cd 7^ 0, decomposes as 

ax + h e , . -ii 1. °-d. , . a 

— h / with e = o and / = — . 

ex -\- d ex + d c c 



This shows that g is output by the composition sequence 
(Mc, Ad, Inv, Me, A/). A more complex example is 



2x _ 1 / 2 

^ (l + a;)2 2 I ^v l+x 



which shows that g is output by the composition sequence 

(Ai, Inv, M_2, Ai, P2, M_i, Ai, M1/2). 
Finally, consider g = log((l + x)/(l — x)). Using 

log ( 1 ■ " ^ 



x-1 

we get the composition sequence (A_i, Inv, M_2, A_2, L). 
Computing the associated power series. Our main 
algorithm requires truncations of the series gi, . . . ,gL asso- 
ciated to a composition sequence. The next lemma discusses 
the cost of their computation. In all complexity estimates, 
the composition sequence o is fixed; hence, our estimates 
hide a dependency in o in their constant factors. 

Lemma 1. If o — (oi,...,ol) is a composition sequence 
that computes power series gi, . . . ,gL, one can compute all 
gi mod x" in time 0(M(n)). 

Proof. All operators in O preserve the precision, except for 
root-taking, since the operator Rfe,ci,r loses r{k~ 1) terms of 
precision. For i < L, define Si — r{k — 1) if Oi has the form 
Rfe,Q,r, Ei = otherwise, and define ul = n and inductively 
rii-i = Ui -\- Ei. Starting the computations with go — x, we 
iteratively compute gi mod a;"' from gi-i mod x"^^^ . 

Inspecting the list of possible cases, one sees that com- 
puting gi always takes time 0(M(ni_i)). For powering, this 
estimate is valid because we disregard the dependency in o: 
otherwise, terms of the form log(fe) would appear. For the 
same reason, 0(M(ni_i)) is in 0(M(n)), as is the total cost, 
obtained by summing over alH. D 

Composition using composition sequences. We now 

study the cost of computing the map Evaln(.,5), assuming 
that g e K[[a;]] is output by a composition sequence o. The 
cost depends on the operations in o. To keep simple expres- 
sions, we distinguish two cases: if o contains no operation E 
or L, we let To(n) — M(n); otherwise, To(n) = M(n) log(n). 

Theorem 1 (Composition). Let o = {oi, ... ,ol) be a 
composition sequence that outputs a series g £ K[[a;]]. Given 
o, one can compute the map Eval„(., p) «n time 0(To(n)). 

Proof. We follow the algorithm of Figure 1. The main func- 
tion first computes the sequence G — gi, . . . , gL modulo x", 
using a subroutine ComputeG(o, n) that follows Lemma 1. 
The cost 0(M(n)) of this operation is in 0(To(n)). Then, 
we call the auxiliary EvaLaux function. 

On input ^,m,n,^,o,G, this latter function computes 
Evalm,n(^, ge). This is done recursively, applying the appro- 
priate associativity rule (Ai) to (A7); the pseudo-code uses 
a C-likc switch construct to find the matching case. Even if 
the initial polynomial A is in K[a::]„, this may not be the case 
for the arguments passed to the next calls; hence the need 
for the extra parameter m. For root-taking, the subroutine 
FindDegrees computes the quantities m; of Eq. (1). 

Since we write the complexity as a function of n, the cost 
analysis is simple: even if several recursive calls are gener- 
ated (fc for feth root-taking), their total number is still 0{1). 



EvaLaux(^, m, n, £, 0, G) 

\f £ = return A mod x" 

£' ^£-1 

switch (o^) 

case(MA): B = Sca\e^,m{A) 

return EvaLaux(i3,m,n,£', 0, G) 
case(A„): B = Shifta,™(A) 

return EvaLaux(i3,m, n, /, 0, G) 
case(Pfc): B = Power™,fc(^) 

return EvaLaux(i3, km — fe + 1, n, £', 0, G) 
case(lnv): B = ReVm(A) 

C = EvaLaux(B, m, n, £' , 0, G) 

return Mul„,„(C,ff]r'" mod a;") 
case(Rfe.a,r): rno, . . . , nik-i — FindDegrees(rn,, k) 

ho^l 

for i = 1, . . . , fc — 1 do 
hi — hhi-i mod x" 

Ao,...,^fc-i = Split^,fe(A) 

for j = 0, . . . , fc — 1 do 

Bi = EvaLaux(Ai, m-i, n, £' ,0, G) 

return Comb„(Bo, • • • , Bk-i, ho,..., ht-i) 
case(E): B = Exp„,„(^) 

return EvaLaux(i3,n, n,i?',o, G) 
case(L): B = Log„, JA) 

return EvaLaux(i3,n, n,^',o, G) 




Eva\(A,n,o) 

G = ComputeG(o, n) 

return EvaLaux(A,n,n, L, 0, G) 





Figure 1: Algorithm Eva I. 

Similarly, the degree of the argument A passed through the 
recursive calls may grow, but only like 0{n). 

Two kinds of operations contribute to the cost: precompu- 
tations of gr^Ti" mod a:" (for Inv) or of l,gi, . . . ,g^'^ mod x" 
for Rk,a,r, and linear operations on A: shifting, scaling, mul- 
tiplication . . . The former take 0(M(n)), since the exponents 
involved are in 0{n). The latter operations take 0(M(n)) if 
no Exp or Log operation is performed, and 0(M(n) log(n)) 
otherwise. This concludes the proof. D 

2.4 Inverse map 

The map Eval„(.,gi) is invertible if and only if g'{Q) 7^ 
(hereafter, g' is the derivative of g). We discuss here the 
computation of the inverse map. 

Theorem 2 (Inverse), x Let o = {oi,...,ol) be a 
composition sequence that outputs g £ K[[a;]] with g'{0) 7^ 0. 
One can compute the map Eva\^^{., g) in time 0{To{ri)). 

Proof. If h is the power series h — '}2i>i hix', with hi^ 7^ 0, 
val(/i) = io is the valuation of h, lc{h) = feip its leading 
coefficient and lt(/i) = higX^" its leading term. We also 
introduce an equivalence relation on power series: p ~ /i if 
g{0) = /i(0) and lt{g ~ g{0)) = lt(ft - /i(0)). The proof of 
the next lemma is immediate by case inspection. 

Lemma 2. For o tn O, if h ^ g and g is in dom(o), then 
h is in dom(o) and o(h) ^ o{g). 

Series tangent to the identity. We prove the proposition 
in two steps. For scries of the form g = x mod x^, it suffices 



to "reverse" step-by-step the computation sequence for g. 
The following lemma is crucial. 

Lemma 3. Let g he in 'Ki\x\], with g = x mod x^ , and let 
o = (oi,...,o_l) be a sequence defined at g. Then o is a 
composition sequence. 

Proof. We have to prove that o is defined at x, i.e., that 
all of oi (a;), 02(01 (a;)), . . . are well-defined. This follows by 
applying the previous lemma inductively. D 

We can now work on the inversion property proper. Let thus 
o = (oi, . . . , o^) be a computation sequence, that computes 
gi, . . . ,gL and outputs g ~ gL, with g = x mod x . We de- 
fine operations 61 , . . . , 6_l through the following table (note 
that we reverse the order of the operations) : 



operation 


Oi 




OL+l-i 


Add 


A„ 




A-a 


Mul 


Ma 




Ml/A 


Powering 


Pfc 


Rk 


.l':(9t-l).val(9i_i) 


Root 


Rfc.Q.r 




Pfc 


Inverse 


Inv 




Inv 


Exp. 


E 




L 


Log. 


L 




E 



Lemma 4. The sequence 6 = (61,..., oz,) is a composi- 
tion sequence and outputs a series g such that g{g) = x. 

Proof. One sees by induction that for all i, 6i_i(- ■ ■ 61(5)) 
is in dom(6i) and 6i{- ■ ■ oi{g)) ~ gL-i- This shows that the 
sequence o is defined at g and that ol{- ■ ■ oi{g)) — x. From 
Lemma 3, we deduce that 6 is defined at x. Letting g be 
the output of 6, the previous equality gives g{g) = x, which 
concludes the proof. D 

Since To = To, and in view of Theorem 1, the next lemma 
concludes the proof of Theorem 2 in the current case. 

Lemma 5. Wtth g and g as above, the map Eval„(.,p) is 
the inverse 0/ Eval„(., 5). 

Proof. Let F be in K[a;]„ and let G = E\/a\„(F,g), so that 
F{g) = G + H, with val(/f) > n. Evaluating at g, we get 
F = G{g) + H{g) = G{g) mod x", since val(p) = 1. D 

General case. Lemma 5 fails when va,\{g) = 0. We can 
however reduce the general case to that where val{g) = 1. 
Let us write g = go + gix -I- • • • , with gi 7^ 0, and define 
g = {g — go)/ gi , so that g = x mod x'^. If o is a composition 
sequence for g, then o = (o, A_j„, Mi/g^) is a composition 
sequence for g, and we have To = To. Thus, by the previous 
point, we can use this composition sequence to compute the 
map Eval^^(.,g) in time 0(To(n)). From the equality 

Eval„(A,g) = Eval„(Scale9i,„(Shiftg„,„(A)),3), 

we deduce 

Evai;:^(^,5) = Shift-go, „(Scalei/gj,„(Evai;;^(A, 5))). 

Since scaling and shifting induce only an extra 0(M(n)) 
arithmetic operations, this finishes the proof of Theorem 2. 

3. CHANGE OF BASIS 

This section applies our results on composition to chanqe 
of basis algorithms, between the monomial basis (a;') and 
various families of polynomials {Pi), with deg(Pi) = i, for 
which we reach quasi-linear complexity. As an intermediate 
step, we present a bivariate evaluation algorithm. 



3.1 Main Theorem 

Let F G K[[a;, i]] be the bivariate power series 

Associated with F, we consider the map 
Eval„(.,F,i) : (ao,...,a„_i) 1-^ I]j<„ ?j(2;)% mod a;". 

The matrix of this map is [Fij]i,j<n. The following the- 
orem shows that for a large class of series F, the opera- 
tion Eval„(., F, t) and its inverse can be performed efficiently. 
The proof relies on a transposition argument, given in §3.3. 

Theorem 3 (Main theorem). Let f,g,h,u,v e ^[[z]] 
be such that 

• g and h are qiven by composition sequences Og and oh; 

• f, u and V can be computed modulo 2" in time T{n); 

• g{0)h{0) = and g'(0),h'{0),u{0),v{0) are non-zero; 

• all coejjictents of f are non-zero. 

Then the series F(a;,t) — u{x)v{t) f(^g{x)h{t)) is well-defi- 
ned. Besides, one can compute the map Eval„(., F,t) and its 
inverse in time 0(T(n) -I- Tog (n) -I- To^(n)). 

Proof. Write / = Efc>o/fe2', 

fl(a;)'° = ^5fe,ia;' and h(t)'' = '^hk,jt^ . 

i>0 j>0 

Since g{0)h{0) — 0, we have that either hkj = for fe > j, 
or gk,i ~ for k > i. Thus, the coefficient F*j of F* is 
well-defined and 



Pi,j = /_, fkgk,ihk,j 



These coefficients are those of a product of three matrices, 
the middle one being diagonal; we deduce the factorization 

Eval„(.,F*,t) = Eval„(.,5)o A„(.,/OoEval^.,ft). 

The assumptions on /, g and h further imply that the map 
Eval„(., F*,i) is invertible, of inverse 

Eva|-i(.,F*,t) = Eva|-*(.,/i)oA„(.,/ri)oEvai;^(.,p). 

By Theorems 1 and 2, as well as Theorem 4 stated below, 
Eval„(., F*,i) and its inverse can thus be evaluated in time 
0(T(n) + Tog(n) -I- To^(n)). Now, from the identity F = 
u{x)v{t)F* , we deduce that 

Eval„(.,F,i) = Mul„,„(.,u)oEval„(.,F*,t)oMul^, „(.,«). 

Our assumptions on u and v make this map invertible, and 

Eva|-'(.,F,i) = Mul^_„(.,fe)oEvai;;'(.,F*,t)oMul„.„(.,a), 

with a{x) = 1/m mod x" and b{t) = 1/v mod i". The extra 
costs induced by the computation of u, v, their inverses, and 
the truncated products fit in 0(T(n) + M(n)). D 

3.2 Change of Basis 

To conclude, we consider polynomials {Pi)i>o in M.[x], 
with deg(Pi) = i, with generating series defined in terms 
of series u,v, f,g,h as in Theorem 3 by 

P = ;^ P{x)e = u{x) v{t) f{g{x)h{t)). 



Eva La ux* (A, m, n, £, o, G) 

\f £ = return A mod x"^ 

£' ^£-1 

switch (of) 

case (Ma): B = EvaLaux*(A, m, n, £', o, G) 

return ScaleA,m(B) 
case (Aa): B = EvaLaux*(A,m, n,/, o, G) 

return Shifta_™(B) 
case (Pfc): B = EvaLaux*(A,mfc - m + l,n,£',o,G) 

return Power^ k{B) 
case (Inv): B = Mul^,„(i,pjr'" mod x") 

C = EvaLaux*(B, m, n, ^', o, G) 

return ReVm(C) 
case (Rfc,c«.r): mo, . . . , ruk-i = FindDegrees(m, fc) 

/io = l 

for i = 1, . . . , fc — 1 do 
hi — hhi-i mod x" 

Ao,...,Ak-i = Comb^A,/io,...,/ifc-i) 

for i = 0, . . . , fc — 1 do 

Bi = EvaLaux*(^i,mi,n,£',o, G) 

return Split^,fc(Bo, • • • , Bfc-i) 
case (E): B = EvaLaux*(yl,n,n,/,o, G) 

return ExpJ„„(_B) 
case (L): B = EvaLaux*(A,n, n, /,o, G) 

return Logl^JB) 




EvalMain'(A, n, o) 

G — ComputeG(o, n) 

return EvaLaux*(A,n,n, L, o, G) 





Figure 2: Algorithm Eval'. 



Corollary 1. Under the above assumptions, one can 
perform the change of basis from {x^)i>o to {Pi)i>o, and 
conversely, in time 0(T(n) + To^ (n) + To^ (^))- 

A surprisingly large amount of classical polynomials fits into 
this framework (see next section). An important special 
case is provided by Sheffer sequences [30, Chap. 2], whose 
exponential generating function has the form 



E 



Hx) 



f = v{t)e' 



•.h(t) 



Examples include the actuarial, Laguerre, Meixner and Pois- 
son-Charlier polynomials, and the Bernoulli polynomials of 
the second kind (see Tables 3 and 4). In this case, if h 
is output by the composition sequence o and v{t) can be 
computed modulo t" in time T(n), one can perform the 
change of basis from (a::*)i>o to {Pi)i>o, and conversely, in 
time 0{T{n) +To(n)). 

3.3 Transposed evaluation 

The following completes the proof of Theorem 3. 

Theorem 4 (Transposition). Leto = (oi, . . . ,ol) be 
a composition sequence that outputs g £ lK[[x]]. Given o, one 
can compute the map Eval^.,^) in time 0(70(71)). 

Proof. This result follows directly from the transposition 
principle. However, we give an explicit construction of the 
transposed map Eval^(.,(;) in Figure 2. Non-linear precom- 
putations are left unchanged. The terminal case £ = 



is dealt with by noting that the transpose of modm,n is 
mod„,m. To conclude, it suffices to give transposed asso- 
ciativity rules for our basic operators. The formal approach 
we use to write our algorithms pays off now, as it makes this 
transposition process automatic. 

Recall that our algorithms deal with polynomials. The 
dual of K[a::]m can be identified with ]K[a;]m itself: to a 
K-linear form £ over ]K[a;]m, one associates '^i^^£{x^)x''. 
Hence, transposed versions of algorithms acting on polyno- 
mials are seen to act on polynomials as well. Remark also 
that diagonal operators are their own transpose. 

Multiplication. In [7], following [19], details of the trans- 
posed versions of plain, Karatsuba and FFT multiplications 
are given, with a cost matching that of the direct product. 
Without relying on such techniques, by writing down the 
rrmltiplication matrix, one sees that Mul^ „(., P) is 

A G K[x]m ^ {AReva+i{P) mod a;"+'') div x"^ G K[x]n, 

if P has degree d. Using standard multiplication algorithms, 
this formulation leads to slower algorithms than those of [7] . 
However, in our usage cases, n, m and d are of the same 
order of magnitude, and only a constant factor is lost. 

Scale. The operator Scale^.n is diagonal; through transpo- 
sition, the associativity rule becomes: 

Eval^,„(^, Mxig)) = ScaleA,m(Eva4,„(^,5r)). (A*i) 

Shift. The transposed map Revjj of the reversal operator 
coincides with Rev„ itself, since this operator is symmetric. 
By transposing the identity for Shift, we deduce 

Shiftl,„(^) = A„(Rev„(Mul^,„(Rev„(A„(A,l/i!)),P)),i!). 

This algorithm for the transpose operation, though not de- 
scribed as such, was already given in [18]. This yields: 

EvalL,„(^,A45)) = Shiftl,^(EvalL,„(A,5)). (A^) 

Powering. The dual map Power^ j, maps A G K[x]^.(^„-i)+i 
to ^o/fc G K[a::]„ (with the notation of §2.1). We deduce: 

Eval^_„(A, Pfc(g)) = Power^^fc(Evalfc(„_i)+i_„(A,g)). (A|) 

Inversion. The transposed version of the rule for Inv is 
Eval^.„(^, Inv(g)) = Rev„(Eval^.„(i\/iul^_„(A,p^""),g)). 

(Ai) 

Root taking. Considering its matrix, one sees that Split^ j. 
maps {Ao, ... , Ak-i) G K[a::]mo x ■ ■ ■ x K[x]m^_i to 

Aoix'') + Ai{x'')x + ■■■+ Ak-i{x'')x''-'' G K[xU. 
Besides, since the map Comb is the direct sum of the maps 

iVlul„,„(.,GO ■.K[x]„ -^K[x]n, 
its transpose Comb^., Go, . . . , Gk-i) sends A G K[a::]„ to 

{Mu\ijA,G,))o<.<k-ienx]n- 
Putting this together gives the transposed associativity rule 

hi = /i' mod a::" for < i < fe 

Ao,. . . , Ak-i = Comb^A, ho,. . . , hk-i) ^ 

Si = Eval^^,„(A,,g) for 0<i < fc ^^^' 

Eval^,„(^, Rfe,a,r(5)) = Split^ .(Bo, • • • , Bfc-i) 



Exponential and Logarithm. From the proof of Propo- 
sition 1, we deduce the transposed map of Exp„ ,j, Log^ „ 
and their associativity rules 

Exp^_„(^) = mod„,™(Shift*_i,„(MultiEval(A„(yl,l/J!)))), 
EvalL,„(A,E(3)) = Exp:„_„(Eval^,„(A,3)); (A^) 

LogL,„(A) = mod„,™(A„(lnterp„(Shiftl,„(A)),i!)), 
Eval^,„(^,L(5)) = Log^,„(Eval^,„(^,3)). (A^) 

4. APPLICATIONS 

Many generating functions of classical families of poly- 
nomials fit into our framework. To obtain conversion algo- 
rithms, it is sufficient to find suitable composition sequences. 
Table 3 lists families of polynomials for which conversions 
can be done in time 0(M(n)) with our method (see e.g. [30, 
3] for more on these classical families). In Table 4, a similar 
list is given, leading to conversions of cost 0(M(n) logn); 
most of these entries are actually Sheffer sequences. Many 
other families can be obtained as special cases (e.g., Gegen- 
bauer, Legendre, Chebyshev, Mittag-Leflter, etc). 

The entry marked by (*) is from [18]; the entries marked 
by (••) are orthogonal polynomials, for which one conver- 
sion (from the orthogonal to the monomial basis) is already 
mentioned with the same complexity in [15, 29]. 

In all cases, the pre-multiplier u{x)v{t) depends on t only 
and can be computed at precision n in time 0{M{n)); all our 
functions / can be expanded at precision n in time 0{n). 
Regarding the functions g(x) and h{t), most entries are easy 
to check; the only explanations needed concern some series 
h{t). Rational functions are covered by the first example of 
§2.3; the second example of §2.3 deals with Jacobi polynomi- 
als and Spread polynomials; the last example of §2.3 shows 
how to handle functions with logarithms. For Fibonacci 
polynomials, the function h{t) = i/(l — t^) satisfies 

l-ft^ 



<-)'^(}^) 



From this, we deduce the sequence for h: 

(Pa, M^i, Ai, Inv, Ma, A_i, Pa, A_i, Ra,a,i, Mi 



/2J 



For Mott polynomials the series h{t) = (1 — \/l — t'^)/t can 
be rewritten as 



h 



l + VT^^t^ 
This yields the composition sequence 

(Pa, M_i, Ai, R2,i,o, Ai, Inv, Ma, A_i, Ra,i,o). 

5. EXPERIMENTS 

We implemented the algorithms for change of basis using 
NTL [32]; the experiments are done for coefficients defined 
modulo a 40 bit prime, using the ZZ_p NTL class (our al- 
gorithms still work for degrees small with respect to the 
characteristic). All timings reported here are obtained on a 
Pentium M, 1.73 Ghz, with 1 GB memory. 

Our implementation follows directly the presentation of 
the former sections. We use the transposed multiplication 
implementation of [7]. The Newton iteration for inverse is 
built-in in NTL; we use the standard Newton iteration for 



square root [10]. Exponentials are computed using the al- 
gorithm of [20]. Powers are computed through exponential 
and logarithm [10], except when the arguments are binomi- 
als, when faster formulas for binomial series are used. For 
evaluation and interpolation at 0, . . . , n — 1, and their trans- 
poses, we use the implementation of [7]. 

We use the Jacobi and Mittag-Leffier orthogonal polyno- 
mials (a special case of Meixner polynomials, with /? = and 
c = —1), with the composition sequences of §2.3. Our algo- 
rithm has cost 0{M{n)) for the former and 0(M(n) log(n)) 
for the latter. We compare this to the naive approach of 
quadratic cost in Figure 3 and 4, respectively. Timings are 
given for the conversion from orthogonal to monomial bases; 
those for the inverse conversion are similar. 




100 200 300 400 500 600 700 800 900 1000 

degree 
Figure 3: Jacobi polynomials. 




100 200 300 400 500 600 700 800 900 1000 
degree 
Figure 4: Mittag-Leffler polynomials. 

Our algorithm performs better than the quadratic one. 
The crossover points lie between 100 and 200; this large 
value is due to the constant hidden in our big-Oh estimates: 
in both cases, there is a contribution of about 20M(n), plus 
an additional M(n)log(n) for Mittag-Leffier. 

6. DISCUSSION 

This article provides a fiexible framework for generating 
new families of conversion algorithms: it suffices to add new 
composition operators to Table 1 and provide the corre- 
sponding associativity rules. Still, several questions need 
further investigation. Several of the composition sequences 
we use are non-trivial: this raises in particular the questions 
of characterizing what functions can be computed by a com- 
position sequence, and of determining such sequences algo- 
rithmically. Besides, the costs of our algorithms are mea- 
sured only in terms of arithmetic operations; the questions 
of numerical stability (for floating-point computations) or of 
coefficient size (when working over Q) require further work. 
Acknowledgments. We thank ANR Gecko, the joint Inria- 
Microsoft Research Lab and NSERC for financial support. 
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/(^) 
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X 


1 - VI - 2t 



Laguerre L° 
Hermite iif„ 
Jacobi Pt'^^ 
Fibonacci F„ 
Euler E^ 
Bernoulli B° 
Mott M„ 
Spread Sn 
Bessel p„ 



polynomial 



2Fi(ii±f±i,H±f±^;/3 + l;^) 



E„>o 7nH4x)t" exp(-t2) exp(^) 

Z^n>0 (/3 + l)„ " y^ 

E„>o ^B-{x)r 
E„>o ^MnWt" 

E„>o S„[x)t 
En>0 ;7T-P"(^)* 

Table 3: Polynomials with conversion in 0(M(n)) 



exp(0) 
exp(«) 
exp(2;) 

cxp(2;) 



generating series 



E„>0 nll^)"* 

E„>o 7nM^)t" 
E„>o^;n&"(^)i"^ 

En>() 771''" (.•^i '^^^ 
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exp(2) 


X 


exp(2) 


X 


cxp(2) 
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h{t) 



Falling factorial (x)„ 
Bell 4>n 

Bernoulli, 2nd kind 6„ 
Poisson-Charlier Cn{x;a) 

(0) 



(*) 



Actuarial a 

Narumi M"^ 

Peters PJ,^'"^ 

Meixner-PoUaczek P„ {x;(l>) (**) 

Meixner m,n,(a;;/3, c) (**) 

Krawtchouk ii'„(a;;p, Af) (**) 

Table 4: 



y 



)K„{x;p,N)r 



1 
1 

t/log(l+i) 

exp(— f) 

exp(/3i) 

riog(l + t)-" 

(l + (l + i)^)-'' 

(l + t^-2tcos0)"^ 

(l + tf 



log(l + t) 
exp(i) — 1 
log(l + t) 
log(l + i/a) 
exp(i) — 1 
log(l + t) 
log(l + t) 



Polynomials with conversion in 0(M(n) log(n)) 



7. REFERENCES 

A. V. Aho, K. Stciglitz, and J. D. Ullinan. Evaluating 
polynomials at fixed sets of points. SIAM J. Comp., 
4(4):533-539, 1975. 

B. K. Alpert and V. Rokhlin. A fast algorithm for the 
evaluation of Legendre expansions. SIAM J. Set. Statist. 
Comp., 12(1):158-179, 1991. 

[3] G. Andrews, R. Askey, and R. Roy. Special functions. 

Cambridge University Press, 1999. 
[4] R. Barrio and J. Peha. Basis conversions among univariate 

polynomial representations. C. R. Math. Acad. Set. Paris, 

339(4):293-298, 2004. 
[5] D. J. Bernstein. Composing power series over a finite ring in 

essentially linear time. J. Symb. Comp.., 26(3):339— 341, 1998. 
[6] D. Bini and V. Y. Pan. Polynomial and matrix computations. 

Vol. 1. Birkhauser Boston Inc., 1994. 
[7] A. Bostan, G. Lecerf, and E. Schost. Tellegen's principle into 

practice. In ISSAC'03, pages 37-44. ACM, 2003. 

A. Bostan, B. Salvy, and E. Schost. Fast algorithms for 

orthogonal polynomials. In preparation. 

A. Bostan and E. Schost. Polynomial evaluation and 

interpolation on special sets of points. J. Complexity, 

21(4):420-446, 2005. 

R. P. Brent. Multiple-precision zero-finding methods and the 

complexity of elementary function evaluation. In Analytic 

Computational Complexity, pages 151-176. Acad. Press, 1975. 
[11] R. P. Brent and H. T. Kung. Fast algorithms for manipulating 

formal power series. J. ACM, 25(4):581-595, 1978. 
[12] P. Biirgisser, M. Clausen, and A. ShokroUahi. Algebraic 

complexity theory, volume 315 of GMW. Springer— Verlag, 1997. 
[13] J. Canny, E. Kaltofen, and Y. Lakshman. Solving systems of 

non-linear polynomial equations faster. In ISSAC'89, pages 

121-128. ACM, 1989. 
[14] D. G. Cantor and E. Kaltofen. On fast multiplication of 

polynomials over arbitrary algebras. Acta Inform., 

28(7):693-701, 1991. 
[15] J. R. Driscoll, J. D. M. Healy, and D. N. Rockmore. Fast 

discrete polynomial transforms with applications to data 

analysis for distance transitive graphs. SIAM J. Comp., 

26(4):1066-1099, 1997. 



[16] M. Frumkin. A fast algorithm for expansion over spherical 

harmonics. Appl. Algebra Engrg. Comm. Comp., 6(6):333-343, 

1995. 
[17] J. g. Gathen and J. Gerhard. Modern computer algebra. 

Cambridge University Press, 1999. 
[18] J. Gerhard. Modular algorithms for polynomial basis 

conversion and greatest factorial factorization. In RWCA '00, 

pages 125-141, 2000. 
[19] G. Hanrot, M. Querela, and P. Zimmermann. The Middle 

Product Algorithm, I. Appl. Algebra Engrg. Comm. Comp., 

14(6):415-438, 2004. 
[20] G. Hanrot and P. Zimmermann. Newton iteration revisited. 

http://www.loria.fr/ zimmerma/papers, 2002. 
[21] G. Heinig. Fast and superfast algorithms for Hankel-like 

matrices related to orthogonal polynomials. In NAA '00, 

volume 1988 of LNCS, pages 361-380. Springer- Verlag, 2001. 
[22] J. g. Hoeven. Relax, but don't be too lazy. J. Symb. Comput., 

34(6):479-542, 2002. 
[23] E. Kaltofen. Challenges of symbolic computation: my favorite 

open problems. J. Symb. Comp., 29(6):891-919, 2000. 
[24] E. Kaltofen and Y. Lakshman. Improved sparse multivariate 

polynomial interpolation algorithms. In ISSAC'88, volume 358 

of LNCS, pages 467-474. Springer Verlag, 1989. 
[25] J. Keiner. Computing with expansions in Gegenbauer 

polynomials. Preprint AMR07/10, U. New South Wales, 2007. 
[26] G. Leibon, D. Rockmore, and G. Chirikjian. A fast Hermite 

transform with applications to protein structure determination. 

In SNC'07, pages 117-124, New York, NY, USA, 2007. ACM. 
[27] Y.-M. Li and X.-Y. Zhang. Basis conversion among Bezier, 

Tchebyshev and Legendre. Comput. Aided Geom. Design, 

15(6):637-642, 1998. 
[28] V. Y. Pan. New fast algorithms for polynomial interpolation 

and evaluation on the Chebyshev node set. Computers and 

Mathematics with Applications, 35(3):125— 129, 1998. 
[29] D. Potts, G. Steidl, and M. Tasche. Fast algorithms for discrete 

polynomial transforms. Math. Comp., 67(224):1577-1590, 1998. 
[30] S. Roman. The umbral calculus. Dover publications, 2005. 
[31] A. Schonhage and V. Strassen. Schnelle Multiplikation grof3er 

Zahlen. Computing, 7:281-292, 1971. 
[32] V. Shoup. A new polynomial factorization algorithm and its 

implementation. J. Symb. Comp., 20(4):363-397, 1995. 



